In [6]:
import sys
!pip install --prefix {sys.prefix} hcp_utils
!pip install --prefix {sys.prefix} nilearn
!pip install --prefix {sys.prefix} nibabel
Looking in indexes: https://pypi.org/simple, https://pypi.ngc.nvidia.com
Requirement already satisfied: hcp_utils in c:\users\adam\anaconda3\lib\site-packages (0.1.0)
Looking in indexes: https://pypi.org/simple, https://pypi.ngc.nvidia.com
Collecting nilearn
  Downloading nilearn-0.9.0-py3-none-any.whl (10.1 MB)
Requirement already satisfied: scipy>=1.2 in c:\users\adam\anaconda3\lib\site-packages (from nilearn) (1.6.2)
Requirement already satisfied: joblib>=0.12 in c:\users\adam\anaconda3\lib\site-packages (from nilearn) (1.0.1)
Requirement already satisfied: scikit-learn>=0.21 in c:\users\adam\anaconda3\lib\site-packages (from nilearn) (0.24.1)
Requirement already satisfied: pandas>=0.24.0 in c:\users\adam\anaconda3\lib\site-packages (from nilearn) (1.2.4)
Requirement already satisfied: requests>=2 in c:\users\adam\anaconda3\lib\site-packages (from nilearn) (2.25.1)
Requirement already satisfied: numpy>=1.16 in c:\users\adam\anaconda3\lib\site-packages (from nilearn) (1.20.1)
Requirement already satisfied: nibabel>=2.5 in c:\users\adam\anaconda3\lib\site-packages (from nilearn) (3.2.2)
Requirement already satisfied: packaging>=14.3 in c:\users\adam\anaconda3\lib\site-packages (from nibabel>=2.5->nilearn) (20.9)
Requirement already satisfied: setuptools in c:\users\adam\anaconda3\lib\site-packages (from nibabel>=2.5->nilearn) (52.0.0.post20210125)
Requirement already satisfied: pyparsing>=2.0.2 in c:\users\adam\anaconda3\lib\site-packages (from packaging>=14.3->nibabel>=2.5->nilearn) (2.4.7)
Requirement already satisfied: pytz>=2017.3 in c:\users\adam\anaconda3\lib\site-packages (from pandas>=0.24.0->nilearn) (2021.3)
Requirement already satisfied: python-dateutil>=2.7.3 in c:\users\adam\anaconda3\lib\site-packages (from pandas>=0.24.0->nilearn) (2.8.1)
Requirement already satisfied: six>=1.5 in c:\users\adam\anaconda3\lib\site-packages (from python-dateutil>=2.7.3->pandas>=0.24.0->nilearn) (1.15.0)
Requirement already satisfied: certifi>=2017.4.17 in c:\users\adam\anaconda3\lib\site-packages (from requests>=2->nilearn) (2021.10.8)
Requirement already satisfied: idna<3,>=2.5 in c:\users\adam\anaconda3\lib\site-packages (from requests>=2->nilearn) (2.10)
Requirement already satisfied: chardet<5,>=3.0.2 in c:\users\adam\anaconda3\lib\site-packages (from requests>=2->nilearn) (4.0.0)
Requirement already satisfied: urllib3<1.27,>=1.21.1 in c:\users\adam\anaconda3\lib\site-packages (from requests>=2->nilearn) (1.26.4)
Requirement already satisfied: threadpoolctl>=2.0.0 in c:\users\adam\anaconda3\lib\site-packages (from scikit-learn>=0.21->nilearn) (2.1.0)
Installing collected packages: nilearn
Successfully installed nilearn-0.9.0
Looking in indexes: https://pypi.org/simple, https://pypi.ngc.nvidia.com
Requirement already satisfied: nibabel in c:\users\adam\anaconda3\lib\site-packages (3.2.2)
Requirement already satisfied: setuptools in c:\users\adam\anaconda3\lib\site-packages (from nibabel) (52.0.0.post20210125)
Requirement already satisfied: numpy>=1.14 in c:\users\adam\anaconda3\lib\site-packages (from nibabel) (1.20.1)
Requirement already satisfied: packaging>=14.3 in c:\users\adam\anaconda3\lib\site-packages (from nibabel) (20.9)
Requirement already satisfied: pyparsing>=2.0.2 in c:\users\adam\anaconda3\lib\site-packages (from packaging>=14.3->nibabel) (2.4.7)
In [10]:
import hcp_utils as hcp
import numpy as np
import nibabel as nb
import nilearn.plotting as plotting
import matplotlib.pyplot as plt
from scipy import signal
import scipy as sci
import random

# load in subject parcel as template cifti
parcelLoc='C:GroCon.dscalar.nii'
parcel=nb.load(parcelLoc)
parcelCort=parcel.dataobj[:,hcp.struct.cortex]

# load in gradient
gradsf='C:hcp.gradients.dscalar.nii'
grads=nb.load(gradsf)
PG=grads.dataobj[0,hcp.struct.cortex]
vox offset (=1.7029e+06) not divisible by 16, not SPM compatible; leaving at current value
vox offset (=1.7029e+06) not divisible by 16, not SPM compatible; leaving at current value
C:\Users\adam\anaconda3\lib\site-packages\nibabel\nifti1.py:590: UserWarning: Extension size is not a multiple of 16 bytes; Assuming size is correct and hoping for the best
  warnings.warn(
vox offset (=1.7029e+06) not divisible by 16, not SPM compatible; leaving at current value
In [2]:
# initialize a percentile membership vector 
percMemb=np.zeros(len(PG.data))
# get percentiles
for p in range(25):
    # every 4 percent bin (25 * 5)
    percThresh_L=np.percentile(PG.data,(p*4))
    percThresh_U=np.percentile(PG.data,((p+1)*4))
    # label cifti verts falling in this perc bin 
    Indices= np.logical_and(PG.data > percThresh_L, PG.data < percThresh_U)
    percMemb[Indices]=p
    
    
    
In [3]:
# plot percentiles/4 to ensure it is working
plotting.view_surf(hcp.mesh.inflated_right, hcp.right_cortex_data(percMemb), 
    bg_map=hcp.mesh.sulc_right)
Out[3]:
In [54]:
# load in template time series to replace data
templateTS=nb.load('C:/Users/adam/template-rest_p2mm_masked.dtseries.nii')
templateData=templateTS.get_fdata()
# extract cortex
templateDataCort=templateData[:,hcp.struct.cortex]

# sampling frequency
TR=.8
fs=(1/TR)

# control for psuedo-high-frequency power from discontinuous (motion-masked) frames
# see https://scipy-cookbook.readthedocs.io/items/ButterworthBandpass.html
from scipy.signal import butter, lfilter
def butter_bandpass(lowcut, highcut, fs, order=2):
    # parameters matched to real-data pipeline
    nyquist = 0.5 * fs
    low = lowcut / nyquist
    high = highcut / nyquist
    b, a = butter(order, [low, high], btype='band')
    return b, a

def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
    b, a = butter_bandpass(lowcut, highcut, fs, order=order)
    y = lfilter(b, a, data)
    return y

# Desired cutoff frequencies (in Hz).
lowcut = 0.008
highcut = 0.09
# create filter
butter_bandpass(lowcut, highcut, fs)

# filter data
templateDataCort = butter_bandpass_filter(templateDataCort, lowcut, highcut, fs)

# extract frequency distribution
# get frequency distribution of all vertices
freqs, density = signal.welch(templateDataCort,fs,axis=0)
# average frequency distribution
avgFreqs=density.mean(axis=1)
plt.plot(freqs,avgFreqs)
plt.xlabel('Hz')
plt.xlim(-0.01,0.25)
plt.show()

# extract cumulative distribution, bin widths (https://stackoverflow.com/questions/13476807/probability-density-function-from-histogram-in-python-to-fit-another-histrogram)
cum_counts = np.cumsum(avgFreqs)
bin_widths = np.repeat(freqs[1] - freqs[0],cum_counts.shape[0])
x = cum_counts*bin_widths
y = freqs[0:]

# interpolate density of frequency distribution
inverse_density_function = sci.interpolate.interp1d(x, y)

# use interpolated density to sample realistic frequencies for simulated waves
randFreqs = np.zeros(100)
for i in range(len( randFreqs )):
    u = random.uniform( x[0], x[-1] )
    randFreqs[i] = inverse_density_function( u )

# plot randomly sampled values
plt.hist(randFreqs,100)
plt.xlabel('Randomly Sampled Values from Distr: Hz')
plt.xlim(-0.01,0.25)
plt.show()
vox offset (=1.70334e+06) not divisible by 16, not SPM compatible; leaving at current value
vox offset (=1.70334e+06) not divisible by 16, not SPM compatible; leaving at current value
C:\Users\adam\anaconda3\lib\site-packages\nibabel\nifti1.py:590: UserWarning: Extension size is not a multiple of 16 bytes; Assuming size is correct and hoping for the best
  warnings.warn(
vox offset (=1.70334e+06) not divisible by 16, not SPM compatible; leaving at current value
In [55]:
### replace template data with sine waves

# initialize a new matrix using template data
SimData=np.zeros(templateDataCort.shape)

# initialize tracker of last-filled TR
LF_TR=0

# for 100 iterations, create a sine wave of realistic frequency
# in each percentile bin
for i in range(100):
    # get realistic frequency
    realisticFreq=randFreqs[i]
    # according to this frequency, how many TRs would this wave take to cycle?
    numTRs=TR/realisticFreq
    # create sine wave of this frequency
    realisticPeriod=np.linspace(-np.pi,np.pi,int(np.round(numTRs)))
    realisticWave=np.sin(realisticPeriod)
    # reverse order of sine wave for EZ interpretation
    realisticWave=realisticWave[::-1]
    # set starting TR for first bin to be 1 TR after where last wave ended
    startTR=LF_TR+1
    # Now distribute it over all gradient bins to match a 20-second procession
    for p in range(25):
        # index vertices in this percentile bin
        thisBin=np.where(percMemb==p)
        # repeat wave to be consistent across all vertices in this bin
        # this is like np.repeat but into a 2d array instead of 1d
        binWave=np.tile([realisticWave],(len(thisBin[0]),1))
        # get ending TR index
        endTR=startTR+len(realisticWave)
        # break loop if end TR is beyond length of time series
        if endTR > templateDataCort.shape[0]:
            break
        # extract frames to index into
        frames=np.arange(startTR,endTR)
        # generate indices that will work to modify SimData
        i1,i2 = np.ix_(frames,thisBin[0])
        # plop realistic wave into next unused TR
        SimData[i1,i2]=np.transpose(binWave)
        # delay 1 TR from previous bin
        startTR=startTR+1
    # record last used TR
    LF_TR=startTR+len(realisticWave)
In [61]:
# need dummy subcortex to save as cifti
OutputData=np.zeros(templateTS.shape)
OutputData[:,hcp.struct.cortex]=SimData
# create cifti from dataframe and templates
timeAxis=nb.cifti2.SeriesAxis(start=0,size=SimData.shape[0],step=1)
spaceAxis=templateTS.header.get_axis(1)
new_img=nb.Cifti2Image(OutputData,(timeAxis,spaceAxis))
filename='C:SimulatedPropagations.dtseries.nii'
new_img.to_filename(filename)
In [60]:
# plot PSD of simulated data
freqs, density = signal.welch(SimData,fs,axis=0)
# average frequency distribution
avgFreqs=density.mean(axis=1)
plt.plot(freqs,avgFreqs)
plt.xlabel('Hz')
plt.xlim(-0.01,0.25)
plt.show()